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Summary 


We proposed to create a single computational code incorporating methods that can model 
both rarefied and continuum flow to enable the efficient simulation of flow about space craft 
and high altitude hypersonic aerospace vehicles. The code was to use a single grid structure 
that permits a smooth transition between the continuum and rarefied portions of the flow. 
Developing an appropriate computational boundary between the two regions represented 
a major challenge. The primary approach chosen involves coupling a four-speed Lattice 
Boltzmann model for the continuum flow with the DSMC method in the rarefied regime. 
We also explored the possibility of using a standard finite difference Navier Stokes solver for 
the continuum flow. With the resulting code we will ultimately investigate three-dimensional 
plume impingement effects, a subject of critical importance to NASA and related to the work 
of Drs. Forrest Lumpkin, Steve Fitzgerald and Jay Le Beau at Johnson Space Center. 

Below is a brief background on the project and a summary of the results as of the end 
of the grant. 
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Calculating Transitional Regime Flows 

Plume flows near spacecraft contain both continuum and rarefied regions. While there 
are many different approaches to modeling each region separately, there is presently no single 
approach that can simultaneously model both regions efficiently and accurately. Having a 
single computational code that can handle both regions is important because the flows may 
be unsteady, there may be multiple interacting regions of rarefied/continuum flow, there can 
be subsonic overlap regions, or the flow computation must reach a steady state through an 
interacting relaxation process during which the different regions exchange information. In 
these cases the use of two separate computer codes is not only awkward and error-prone but 
also inefficient due to the need to swap and convert data between the codes after every time 
step. This is particularly so because the codes may not represent the fluid in the same way. 

In a fully continuum flow the integration of the Navier Stokes (NS) equations is a well 
established practice [1]. For lower density flows, however, the NS equations no longer apply. 
The Burnett equations may be used to extend the applicability of the NS equations further 
into the rarefied regime [2]. However, the most practical approach for modeling rarefied 
flow is the Direct Simulation Monte Carlo (DSMC) method of Bird [3] in which Lagrangian 
molecules are collided and moved to model the flowing gas. The computational cost in 
DSMC is proportional to the number of molecules simulated. To model transitional flow, 
for example continuum-to-rarefied plume flow expansions, the DSMC has been used all the 
way into the nozzle throat. In fact, this is one approach studied at JSC by Jay Le Beau 
in the Aeroscience Branch of the Navigation, Control and Aeronautics Division. However, 
except for very small plumes [4], it has been impractical to use full resolution for the DSMC 
simulations inside the nozzle. Thus, for purposes of computing the far field low density flow, 
the continuum flow was only crudely modeled by DSMC. 

Merging Continuum and DSMC Codes 

We are developing a hybrid computational code which can handle unsteady transitional 
flows containing both fully rarefied and fully continuum regions (i.e., the flow field in and 
around a thruster during startup). To do this we are merging the DSMC and a continuum 
procedure into a single routine. This work has centered around developing the interface 
between the two different representations of a flow: macroscopic parameters (U, p , T etc.) 
in the continuum flow and the microscopic particle velocity distribution functions in the 
rarefied flow. For unsteady and subsonic flows it is essential that the flow information be 
free to propagate in either direction, i.e., rarefied to continuum or continuum to rarefied. 
It is also desirable that both types of code be compatible with future vector and parallel 
supercomputer architectures. 

We have been working on the project with PhD students who have passed the doctoral 
qualifying exam and have had JSC funding since June, 1995. First, we explored a range of 
code combinations since it was not obvious which types of codes should be merged. One com- 
bination appeared well suited to the problems of interest: the Lattice Boltzmann Equation 
(LBE) discrete-velocity approach [5,6] for use in the continuum region and the DSMC ap- 
proach in the transitional and rarefied regions. The advantages of using the discrete-velocity 
gas model are that then the representation of the flow in terms of velocity distribution func- 
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tions and particles is similar in both flow regimes. In addition, it may be possible to compute 
flows more quickly using less memory than with an ordinary finite difference representation 
of the continuum flow. 

Recently, the continuum Euler equations have been solved using a discrete- velocity distri- 
bution function model based on a local thermodynamic equilibrium assumption. This LBE 
method is termed the Adaptive Discrete Velocity (ADV) method of Nadiga [7]. The approach 
uses a flux-splitting/TVD scheme based on Nadiga and Puffin’s Flux-splitting Equilibrium 
Method [8]. The ADV approach obviates the need for the expensive computation of the 
Lattice Boltzmann collision integral by assuming an equilibrium distribution function based 
on macroscopic properties at each cell center. However, the method is also unique in that 
it retains a non-equilibrium distribution function at cell boundaries because of the superpo- 
sition of fluxes arising from neighboring cells with possibly different properties. The ADV 
approach appears to allow an efficient coupling between continuum and rarefied regions by 
combining a quasi-particle approach with a true particle simulation. 

Hence, we have chosen to merge the DSMC and ADV/LBE approaches because: a) 
DSMC is the most common, flexible, and well tested solver for the rarefied flow regime, 
b) the ADV approach provides a suitable representation of discrete distribution functions 
to interface with the DSMC, c) the work greatly benefits from on-going collaboration with 
Dr. Nadiga at the Los Alamos National Laboratory, d) the ADV approach is a reasonably 
efficient and stable approach to solving the continuum flow, and e) both methods maintain 
similar references to distribution functions and proceed through related move-collide-move 
and flux-equilibrate-flux sequences. On the other hand, finite difference/ volume NS solvers 
are now industrial standard approaches having well documented features and limitations. 
Thus, we are developing an interface that is not entirely specific to the ADV and DSMC 
codes and should be more broadly applicable to merge other NS solvers with DSMC. 
Results 

As part of this grant we have developed suitable versions of the DSMC and ADV codes 
and have created a reasonably satisfactory but simple interface. The interface consists of 
a layered set of reservoir cells, buffer cells and ghost cells (Fig. 1). The particles which 
enter the DSMC region are generated from the equilibrium ADV properties in the two 
ADV cells bounding the DSMC region. The boundary cells in the ADV domain (ADV 
buffer cells) similarly utilize the macroscopic properties in the two bounding DSMC cells 
to obtain the DSMC-to-ADV fluxes. Because the ADV method is second order in time 
and uses half-step properties, a form of half-step DSMC is done in the boundary cells. To 
obtain smooth time-accurate inputs to the ADV region, there is a portion of the DSMC 
region along the boundary in which ‘ghost-cells’ are used. Within the ghost-cells essentially 
an n— parallel DSMC computation occurs in which the n ghost levels compute statistically 
somewhat independent problems. 1 Particles entering the ghost region are cloned n times 
while those departing for the ordinary DSMC side are destroyed with probability n/(n + 1). 
Those which leave for the ADV side are lost altogether. Hence, we may compute with 40 

x The problems are not truly independent because each ghost level has cloned particles as inputs. 
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particles per cell in most of a DSMC region while, with 10 levels of ghost-cells (n = 10), 
the ADV buffer cells see a random noise level closer to that found with 400 particles per 
cell. Presently we use a rectilinear Cartesian grid in both regions although the cell sizes may 
differ across the boundary. A set of grid maps are used to delineate the different regions and 
cell types. 

Consider some example results illustrating the use of the hybrid code. See the figure 
caption page for details of each simulation. Figure 2 shows the unsteady passage of a normal 
shock wave through the interface. The passage occurs smoothly with no apparent reflected 
waves and the proper shock jump conditions are maintained. Notice the change is shock 
slope between the two methods: whereas DSMC fully resolves the wave, the Euler solution 
of the ADV method produces a wave whose thickness (roughly 3 cells) is controlled by 
second order numerical viscosity. Passage of rarefactions is equally uneventful. Figure 3 
illustrates the passage of a shear “wave” through the interface. The problem solved is the 
strongly compressible Rayleigh problem: a diffusely reflective flat plate was abruptly moved 
in its own plane. This generated weak shocks and rarefactions which quickly left the scene. 
Much later, as seen in the figure, the diffusively growing shear layer gradually penetrates 
the interface and experiences a slope discontinuity due to the inconsistent nature of the 
problem (nearly an Euler/inviscid solution on the ADV side). This example illustrates that 
the interface should remain within the nearly inviscid flow to avoid such errors. Figure 4 
illustrates the use of DSMC to capture a rarefied flow effect within an otherwise continuum 
flow. A Mach 3 shock wave is captured by DSMC while ADV computes the upstream and 
downstream regions. The quieting effect of the ghost cells is evident. Figure 5 presents 
preliminary results for a simple plume impingement problem. A small plume, generated 
by a 10 to 1 pressure difference across a 2D narrow slot, strikes a flat plate 6 slot heights 
downstream. This problem demonstrates several features of our approach: (1) solid surfaces 
on which boundary layers may be important are within the DSMC domains, (2) multiple 
detached domains can be used, (3) interfaces can be placed in sub- or supersonic flow, (4) 
oblique wave passage through the interfaces is uneventful, (5) continuum code can be used 
to handle subsonic domain boundaries, and (6) reasonable solutions can be obtained with 
minimal ensemble averaging. A video tape containing the unsteady development of this flow 
has been sent to Dr. Lumpkin. 

Finally, we point out that a detailed description of the work was presented at the 35th 
Aerospace Sciences meeting in Reno, NV in January, 1997 as paper 97-1006. This paper is 
being prepared for submission to a journal and has been sent to Mr. Cones by surface mail. 
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Figure Captions 

Figure 1. Interface Configuration. 

The follwing information applies to all four hybrid ADV /DSMC simulations in the follow- 
ing figures. “Test conditions” indicate atmospheric density at 80 km altitude (9.2e-6 kg/m 3 ) 
and a temperature of 500 K. The grids employ equi-sized ADV and DSMC cells measuring 
0.8 rest mean free paths per side for figures 2,4,5 and 0.5 rest mean free paths per side for 
figures 3. DSMC hard sphere particles have the mass of nitrogen molecules. All solutions 
are advanced with time steps measuring 1 x 10“ 6 seconds and all were done on a DEC Alpha 
3000/700 workstation. 

Figure 2. Moving Normal Shock Wave. Figure 2 shows 11 realizations from a hybrid 
ADV/DSMC simulation of a normal shock wave passing through the interface. 

The whole domain consists of 80 DSMC cells adjacent to a specular wall on the right 
hand side of the domain, and 52 ADV cells to the left of the interface. We further subdivide 
the DSMC region into 20 ghost cells, n = 5, near the ADV/DSMC interface. Each cell has 
200 simulated molecules upstream of the shock. Thirty minutes were required to advance 
the computation through the 1000 time steps shown. 

Flow initially moves uniformly toward the solid wall at 400 m/s (Mach 1.6). We initialize 
the constant upstream inflow to test conditions. Periodic boundary conditions are applied 
in the y direction. 

The shock thickness spans about 12 rest mean free paths in the DSMC region, whereas 
in the ADV region it is thinned considerably and has a finite thickness only due to the effect 
of numerical viscosity. The solution yields the correct Rankine- Hugoniot density ratio. In 
the ghost cells the five-fold increase in ghost particles produces a clear reduction in noise. 

Figure 3. Compressible Rayleigh Problem. Figure 3 shows tangential velocity profiles 
beside a flat plate suddenly moved in its own plane. 

This flow demonstrates the nature of shear waves passing through the ADV /DSMC 
interface. We compare two hybrid DSMC/ADV solutions with differing DSMC domains. 
The whole domain consists of 400 cells. One computation (dashed curve) uses 280 DSMC 
cells adjacent to the specular wall whereas the second computation uses only 80 DSMC 
cells. Both DSMC domains use 20 ghost cells, n = 5, near the ADV/DSMC interface. Each 
DSMC cell is initialized with 200 simulated molecules. The solution is advanced through 
10000 iterations, three of which are plotted after 200, 6000 and 10000 steps. Thirty minutes 
of computation time were needed. 

Periodic boundary conditions are applied in the vertical direction while in the horizontal 
direction an absorbing boundary condition is imposed at the left side of the ADV domain. 
The plate at the right side of the DSMC domain is suddenly moved in its own plane to 2000 
m/s (Mach 4) at t = 0. 
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There exists a small velocity slip adjacent to the wall. The full DSMC simulation captures 
the type of decay to the farfield value (zero u velocity) similar to that in the incompressible 
Rayleigh problem. The choice of collision cross section determines the viscosity in the DSMC 
method whereas only numerical viscosity affects the ADV result. The discontinuity in vis- 
cosity at the interface leads to thinning of the viscous region and a corresponding increase 
of the velocity gradient in the ADV. 

Figure 4- Steady Shock. Figure 4 shows 47 overlaid realizations from a simulation of a 
stationary, normal shock wave. 

We initialize each of the 62 DSMC cells with 800 simulated molecules upstream and 
2400 downstream of the shock. Two DSMC border ghost zones (12 cells wide, n = 5 levels 
per cell) provide smooth data to the surrounding ADV regions. The solution is advanced 
through 5000 time steps and took 30 minutes. 

Rankine-Hugoniot equations determine the boundary conditions in the x direction while 
periodic boundary conditions are applied in the y-direction. We initialize the Mach 3 up- 
stream inflow conditions to test conditions. The simulation yields the correct perfect gas 
density ratio (3.0), as expected and the correct shock shape. The shock half-slope thickness 
spans approximately 4 rest mean free paths. 

The effect of ghost cells is clearly visible at either end of the DSMC domain. The five- 
fold increase in ghost particles produces the expected 56% noise reduction. The exponential 
behavior of the noise in the ghost zone near the left interface arises due to cloning of incoming 
equilibrium particles through all n levels. Collisions in the cells downstream of the interface 
destroy the cloned identity, leading to noise reduction in the averaging procedure. Since the 
upstream flow is highly supersonic, noise propagation occurs only in the subsonic downstream 
ADV region. 

Figure 5. Supersonic Jet Flow Impinging on a Vertical Flat Plate. Figure 5 shows 
pressure contours (Fig. 5a) from a 2D calculation of supersonic jet flow impinging on a 
vertical flat plate. 

Two separate DSMC patches, embedded in the ADV domain, contain the slit vertical 
walls and the target plate, respectively (Fig. 5b). The whole domain (152*84 cells) consists 
of 2706 DSMC cells and 10062 ADV cells. We model the slit with two vertical specularly 
reflecting walls on either side of a centered opening 9.6 rest mean free paths (12 cells) high. 
The specularly reflecting target plate measures 4.8 rest mean free paths (6 cells) high and is 
positioned 60 rest mean free paths (75 cells) to the left of the slit. 

To create a starting pressure ratio of 10:1, we initialize the DMSC cells to the left of the 
slit with 7 simulated molecules and those to the right with 70 molecules. The flow develops 
from the intial state of rest and test conditions assigned to the low pressure region. The 
simulation has a Reynolds number of 10 based on the upstream speed of sound and the plate 
height. 

Periodic boundary conditions are applied in the vertical direction whereas absorbing 
boundary conditions are enforced in the horizontal direction. The solution is advanced to 
steady flow through 20 realizations of 100 time steps each. Fig. 5a shows pressure contours 
averaged over the last 4 independent realizations. The total run time was 2 hours and 15 
minutes. 

Choked flow exists at the throat. Pressure contours show the formation of an expansion 
region bounded in the vertical direction by shear layers (free boundaries) and terminated by 
a diffuse normal shock at x=0.24 meters. Subsequent expansion results in supersonic flow 
which causes the bow shock ahead of the vertical plate. Two recirculating regions arise on 
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either side of the expansion region. 
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